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We introduce and study a da Vinci Fluid, a fluid whose dissipation is dominated by solid friction. 
We analyse the flow rheology of a discrete model and then coarse-grain it to the continuum. We find 
that the model gives rise to behaviour that is characteristic of dense granular fluids. In particular, 
it leads to plug flow. We analyse the nucleation mechanism of plugs and their development. We 
find that plug boundaries generically expand and we calculate the growth rate of plug regions. In 
systems whose internal effective dynamic and static friction coefficients are relatively uniform we 
find that the linear size of plug regions grows as (time)^''^. The suitability of the model to granular 
materials is discussed. 

PACS numbers: 47.57.Gc, 45.70.Mg, 62.20.Qp 



I. INTRODUCTION 

The flow of non-cohesive granular matter has focused 
much attention due to both a broad technological rel- 
evance and deep scientific challenges. The importance of 
modeling granular flow cannot be overemphasized - par- 
ticulate transport is important to dry chemicals, phar- 
maceutical granules, agricultural grains, cereals, pebble 
beds in nuclear reactor and more. Of particular interest 
is flow in hoppers and silos, commonly used for storage 
and discharge. Dense particulate flows are also relevant 
to many natural phenomena: avalanches of rocks and 
snow, cratering, quicksand dynamics and dune locomo- 
tion. 

Flow of dilute systems is straightforward to model, us- 
ing pair inelastic collision theories In many cases, 
however, the flow is dense and even the concept of col- 
lision is not well-defined when grains are in prolonged 
contact. Many attempts to model dense flow p^-fT^I fol- 
low a trend which favors increasing model complexity to 
accommodate addition of physical mechanisms. While 
this helps in fltting experimental measurements, it makes 
models less general and difficult to extract generic insight 
from. 

Here we focus on a minimal model - a fluid governed 
solely by solid friction. Solid friction plays a particularly 
important role in the fiow of dense non-cohesive partic- 
ulates - it provides a mechanism to dissipate mechanical 
energy as grains rub against neighbours, and to store it 
in intra-granular degrees of freedom. We call this a da 
Vinci fiuid (dVF), after da Vinci's pioneering work on 
solid friction [isj , which preceded the more known works 
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of Amontons [T3l and Coulomb [11]. In a dVF, volume 
elements follow Newton's second law, but they interact 
via the da Vinci - Amontons - Coulomb law of solid fric- 
tion. We show that this model captures correctly one 
of the most important aspects of flow of dense granular 
materials - the formation and growth of plug flows. We 
present exact solutions under simplified conditions that 
provide insight into the effect of this mechanism on the 
macroscopic rheology. 

In the following we first formulate the fiow equations for 
laminar flow in a model discrete system. We then solve 
the equations for several simple cases. Next, we use the 
gained insight to obtain continuum equations. Finally, 
we derive the dynamics of growth of simple plug regions. 
We construct our model for laminar flow of dense gran- 
ular material by dividing the system into layers, each 
of many grains, perpendicular to the direction of flow. 
Individual grains across the common interface between 
adjacent layers interact via normal and friction contact 
forces. These give rise to mean normal and tangential 
friction forces between the layers. Following experimen- 
tal and numerical observations fiH|[i3|' '^^ assume that 
the threshold for relative motion between layers occurs 
once the ratio of the friction force to the normal force 
exceeds an effective static friction coefficient /Ug. Once 
relative motion has been established, the friction force is 
velocity-independent [isll and is proportional to the nor- 
mal force with a dynamic friction coefficient /i^ < /^s • 
The discrete system is shown in figure 1. It comprises 
a set of N parallel layers of width d, all moving in the 
z-direction. The layers are confined between two bound- 
ary plates and oriented in the direction of gravity. A 
pressure P is applied to the boundary plates, as shown 
in figure 1. Effective dynamic and static friction coeffi- 
cients are assumed both between the layers and between 
the outermost layers and the boundary plates. The focus 
on planar layers is mainly for convenience - the following 
analysis is general and applies to any uni-directional fiow, 
such as parallel streamlines in a pipe. For simplicity, we 
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assume that the layers have a uniform mass density p, 
but the model can be extended to non-uniform density, 
in which case the local effective friction coefficient due to 
friction between grains in adjacent layers, depends on the 
local density and hence must be position dependent. In 
the present article we consider, however, only cases where 
the initial density is uniform and remains so throughout 
the dynamics. More general cases, where the friction co- 
efficient depends on the inertial number / [l3| [HI [HI , 
will be addressed in a later report. 

When moving relative to the boundaries, the layers mo- 
tion in the z-direction is given by 



pdvn = pdg + < 



Pn+i,n - ^J■dP if n==l ; 

if l<n<N 

Pn~l,n - Md-P */ n=N , 



(1) 



where n — 1,2, N indexes the layers and Pn+i,n is the 
friction force per unit area that layer n-t-1 applies on layer 
n. Note that the above equations are correct whenever 
g > HsP/ipd), regardless of the relative velocity between 
the boundary and the outer layers. By Newton's third 
law, this force is equal and opposite to the force that 
layer n applies on layer n + 1, Pn+i^n = -Pn,n+i- 
The velocity of the centre of mass of all the layers is 
vcM = X^n'*'"/^ and, by summing the equations, it is 
straightforward to verify that the acceleration of the cen- 
tre of mass is 



vcM = g 



2pdP 
pL 



(2) 



where L is the distance between the stationary boundary 
plates. The friction forces on the right hand side of eqs. 
([T]) follow the da Vinci - Amontons - Coulomb law and 
depend non-analytically on the relative velocity between 
neighboring layers. 

Noting the significance of plug flow, we wish to determine 
the consistency of such solutions with the threshold na- 
ture of the friction forces. We therefore consider first a 
case where all the initial velocities and accelerations of 
the layers are equal, given by eq. 0. Calculating the 
friction force per unit area between layers n and n + 1, 
we find 



Pn+l,n = P-dP 1 



2n 

'n 



(3) 



decreasing linearly with distance from the left boundary. 
From ([3]), all the inter- layer friction forces are below the 
threshold value PsP, i-C- the layers cannot slide relative 
to one another and they move together. Thus, there ex- 
ists a uniform plug flow solution where all the layers fall 
as one rigid body. 

Next, let every layer move initially relative to both its 
neighbors, with an overall velocity profllc that increases 
from left to right. 



(4) 



Since layer 2 falls faster it applies on layer 1 a friction 
force (per unit area) p2,i — PdP in the z-direction. This 
force is cancelled by the friction force on 1 from the left 
boundary and 1 accelerates at ai = g. Similarly, the 
friction forces on every layer n < N cancel out, leading 
to all these layers accelerating initially at g. Thus, for a 
short time t, before any threshold is reached, the velocity 
profile for 1 < n < remains unchanged, Vn+i — Vn = 
constant. 

This solution, however, is unstable. Layer N experiences 
decelerating friction forces from both layer N — 1 and 
the right boundary (figure 2) and its acceleration is g — 
2pdP/pd < g. Therefore, layer — 1 gains on layer A^ 
and, after a time 



2pP 
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gTi. At that 



their velocities match at v(ti) = 
moment pn,n-i vanishes and the two layers move as one, 
forming a nuclear plug. From ri on, eqs. ([Ij need to be 
solved afresh with a new velocity profile 



if n<N ; 
if n=N . 



(6) 



The acceleration of the nuclear plug is g — p^P/ pd < g. 
Layer N — 2 then starts gaining on the plug, catching 
up with it at T2 = pd{v'^_^ — w^^j)//^^^- At that mo- 
ment layer N — 2 joins the plug and the three layers move 
as one. This involves instantaneous readjustment of the 
friction forces between the layers within the plug: pn,n-i 
becomes —pdP/5 and pn-i,n-2 drops below the thresh- 
old to pd.P/'i- The process of velocity matching and read- 
justment of inter-layer friction forces in the growing plug 
region (PR) continues from right to left until it reaches 
layer 1. From then on, the entire system is a plug, mov- 
ing as one rigid body. Thus, the flow converges to the 
flrst case discussed above. 



The time interval, r,„, 
N — m and N — m 



between the moments that layers 
1 join the plug is 



pd{vl 



,,0 

-'N-m+ 



l) 



2pdP 
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and the acceleration of the right m-layer plug during this 
interval is g — 2pdP/'nipd. The time it takes the system 
to converge to a global plug is 



N 



(8) 



and then it accelerates at 5 — 2pdP/pL. The evolution 
of the velocity profile is shown in figure 1. 



FIG. 1: The evolution of a graded velocity profile (dashed 
line), shown at successive times from top to bottom. A PR 
forms at the right boundary and expands into the fluid. 



Next, consider an initial arbitrary non-monotonic veloc- 
ity profile w^, containing no PR. Layers whose velocities 
are neither a local maximum nor a local minimum of 
move faster than one neighbor and slower than the other. 
Hence, the friction forces on them cancel out and they 
accelerate at g. When the velocity of the layer is a lo- 
cal maximum of the velocity profile (figure 2), both its 
neighbors slow it down and it accelerates at 5 — 2ij,aP. 
Similarly, when a layer's velocity is a local minimum, its 
neighbors act to accelerate it at g + 2^dP- Both these 
cases lead to plug nucleation. Thus, the catch-up dy- 
namics described above generate nuclei of plugs at all 
the extrema of wJJ and these plugs subsequently expand. 
When two expanding PR's meet they merge and continue 
as one plug. The merging dynamics is interesting in its 
own right, but it is downstream from this discussion (but 
see comment below) . The expansion and merging of PR's 
continues until the entire system is one plug. Thiis, the 
global plug flow discussed above is a stable fixed-point 
solution of a whole family of laminar flows. 
The coalescence of PR's is far from simple. When the 
boundaries of two PR's collide they are at different ve- 
locities and accelerations and the velocity matching be- 
tween them is accompanied by readjustments of inter- 
layer forces within each PR. This gives rise to interest- 
ing dynamics that depend on intra-plug stress relaxation 
rates, which are not described by our equations. How- 
ever, by assuming that the readjustment of forces is in- 



FIG. 2: The evolution of a velocity profile containing a local 
maximum. The velocity profile (dashed line) is shown at suc- 
cessive times from top to bottom. The fastest layer is slowed 
down by friction with both its neighbors, who eventually catch 
up with it. A plug then nucleates and expands outwards of 
the maximum point. 



stantaneous, we can compute the intra-plug inter-layer 
forces and work out explicitly the catch-up dynamics all 

the way to global plug fiow. 

The above convergence to a uniform plug state is not 
necessarily the rule. For example, consider a flow in a 
pipe, tilted at an angle a to the horizon (figure 3). The 
fluid is supported by a bottom plate, whose friction with 
the fluid is the same as the internal friction, and is con- 
flned by a top frictionless plate. Signiflcantly, now the 
inter-layer pressure varies with position, since upper lay- 
ers weigh down on lower ones, leading to varying inter- 
layer slippage thresholds. 

It is straightforward to show that, if all layers are initially 
at rest imdcr a uniform boundary pressure Pq, n top 
layers slide down as one, while the bottom N — n layers 
remain at rest. The value of n depends on the angle a 
and on the friction coefficient /Xg, it is the smallest integer 
that satisfies the relation 



g pd{cos a — IJ.S sin a) 

The n-layer PR does not expand to span the entire sys- 
tem because, once it starts sliding, the friction coeffi- 
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FIG. 3: A system of A'^ identical layers tilted to the direction 
of gravity g. 



cient between it and the n + 1th stationary layer drops 
to fid, reducing the tangential force. This demonstrates 
(demonstrating) that the dVF equations accommodate 
stable flow solutions of localized PR's. The analysis can 
be extended to arbitrary initial velocity profiles. In par- 
ticular, it is possible to identify the thresholds for inter- 
layer relative motion and the times that layers join the 
PR. These give rise to rich patterns of sliding regions, 
which will be explored elsewhere. 



We now use the insight from discrete systems to formu- 
late continuum flow equations. We first rewrite the dis- 
crete equations to take into account the non-differentiable 
form of the friction forces. Then we make the equations 
continuous by tending the layers' widths to zero, mak- 
ing them into streamlines. The dynamics depend non- 
analytically on velocity differences between layers and 
the equation of motion for layer n can be written as 



1 

pd 
pd 



Vn+l - V„ 



P„+iSign 



Vn+l - Vn 



Vn - Vn-1 
OL I ; I Vn,n-\ 



P„_iSign 



(10) 



where we define, respectively, Sign(u) = —1,0,1 when 
M<0,M = 0,M>0 and a{u) = 0, 1 when u ^ 0, u = 0. 
The expression in the first brackets on the right hand side 
describes sub-threshold friction forces when layers move 
together. The second brackets contain the contribution 
from friction forces between neighbour layers at different 
velocities. Sub-threshold friction forces satisfy p„+i,„ -f 
Pn-i,n = FpB., where FpR is the total force on the PR 



boundary. 

To extend eqs. (jlOp to the continuum, we define a contin- 
uous coordinate, x = nd, running normal to the bound- 
aries from left to right (figure 1) and take the limits d ^ 
and TV — > cx), such that Nd = L. For generality, we let 
P = P{x) be a function of position. In this limit we 
obtain 



J 



dv 
dt 



ld_ 

p dx 



p{x) a 



dv 
dx 



p dx 



P{x) Sign 



dv 
dx 



(11) 



r 



where p{x) is the friction force per unit area between where 9f /9a; changes discontinuously from zero to a finite 
streamlines. As in discrete systems, this flow is unstable value, 
to nucleation and growth of PR's, whose boundaries are 

We next obtain an equation for the boundary growth. 
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Consider a fluid, containing at least one PR, flowing in 
the z-direction under gravity g with a velocity profile 
Vp{x) at t = 0. Consider a PR of width w between xi 
and Xr — xi + w moving at velocity uq- The left and 
right streamlines just outside the PR move, say, slower 
than the plug and, for concreteness, we choose their ini- 
tial velocities as vi < Vr < uq. The streamlines apply a 
deccelerating force of on the PR and, after a short 

time t, the velocities of the PR and its boundary stream- 
lines are 



Vrit) 

u{t) 

Vlit) 



vr{0) +gt 



Uq 
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pw 



(12) 



vi{0)+gt 



The right streamline catches up with the plug first after 
Tr = pw [uq — Vr{0)] /2/idP. Dividing both sides of this 
expression by d and taking the continuum limit d — ?■ 0, 
the expansion of the PR to the right follows 



dt 



pw 



dv 
dx 



(13) 



where (dv/dx)^ is the gradient of the z-directed velocity 
at the right boundary. Combined with a similar expan- 
sion to the left, we find that the plug broadens at a rate 



dw 
'dt 



2pdP 
pw 



dv 
dx 



dv 
dx 



(14) 



Suppose the initial flow velocity profile is analytic with 
a local maximum at a; = 0, v{x, t = 0) = uq — Cx^ + . . .. 
This gives birth to a PR that broadens, according to (ITil) . 
as 



dw 



pCw^ 



2pdP _ 
■ipC ' 



1/3 



(15) 



This broadening as t^^^ is generic, occurring in different 
geometries and higher dimensions, as will be reported 
elsewhere |20||. 



II. CONCLUSION AND DISCUSSION 

To conclude, we have presented a minimal description 
of flow of dense granular matter, modeled as a da Vinci 
fluid. The dVF is governed by Newton's equations with 
normal contact forces and drag due to solid friction be- 
tween volume elements. Consequently, accelerations can- 
not be ignored, as in many existing models (e.g. [2^ . We 
first analysed the behavior of discrete systems, showing 



that the threshold nature of the da Vinci - Amontons - 
Coulomb friction law gives rise to rich dynamics. Most 
notably, we have shown that the flow is unstable to nucle- 
ation of plug regions, wherein the fluid flows at a uniform 
(position-independent) velocity. We have found that, un- 
der some conditions, once a PR has nucleated, it may 
expand. We have identified a growth mechanism, which 
we call catch-up dynamics, and illustrated the dynamics 
explicitly in several examples that provide clear physical 
insight. 

We have then extended the formulation to the contin- 
uum, where the discrete layers become streamlines whose 
drag is governed by the da Vinci - Amontons - Coulomb 
friction. We find that, under appropriate conditions PR's 
expand, in which case their linear size grows proportion- 
ally to i^/'^. In these cases PR's expand and coalesce 
until the entire fluid flows as one plug. We have also 
shown that, under some boundary loading, PR's may re- 
main flnite. It is interesting that dense clusters forming 
in granular gases also grow as t^^^ plj . 
The normal hysteresis associated with the different values 
of ps and pd, hardly plays a role in the model. This is 
because the value of p^ is only implicit as a threshold 
in the PijS in eq. pop . The hysteresis would come into 
play when intra-plug layers start moving relative to one 
another, but there is no mechanism to initiate internal 
motion. Forces are applied to plugs only through their 
boundaries, via dynamic friction. 

The model has several advantages. First, it is minimal, 
involving only inertia and solid friction, and hence it is 
straightforward to interpret physically. In particular, the 
model alleviates the need to resort to an energy equation 
since the energy, for this type of fluid, can be computed 
directly from the equations. Indeed, the irrelevance of 
energy as a variable in linear Couette flow of granu- 
lar materials has been discussed recently by Kumaran 
[1^ . Second, it gives rise generically to plug flow, a phe- 
nomenon often observed in flow of dense granular fluids, 
via a tractable physical mechanism. The simplicity of the 
model makes it useful for gaining insight into the forma- 
tion and development of plug regions, and in particular 
into their peculiar rate of expansion. 
It should be noted that, in contrast to many existing 
models, which assume a steady state flow of constant ve- 
locity profile, there is no such assumption in the dVF 
model. For example, under the pure gravitational flow 
shown in figures 1 and 2, the steady state consists of a 
plug flow moving at a constant acceleration. This is a di- 
rect consequence of considering only solid friction drag. 
This also sets the range of validity of the model because, 
as local velocities increase under the acceleration, the 
mean number of contacts per grain at any moment de- 
creases. Once this value drops below one, the material 
can be considered dense no longer and the flow crosses 
over to what is known as the inertial regime. Thus, this 
model, while applying strictly to dense flows, also shows 
how the fluid can exit the dense regime. 
To keep the model minimal, we have assumed uniform 
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friction coefScients throughout the materiaL This simpU- 
fication gives good insight into the phenomenon of plug 
formation and growth, but it prevents capturing some 
phenomena. Notably, it does not predict shear banding 
and boundary pressures that vary with boundary shear 
rates (e.g. as in gravitational chute flows, where, given an 
inertial number / j"l8| . the pressure varies as the square 
of the shear rate). It also leads directly to convergence 
of the flow to a plug that spans the entire volume. Re- 
alistically, the local effective friction coefficient may be 
a function of the local density, as well as of microscopic 
grain-scale details. For example, the local value of /i can 
be suppressed considerably by rolling of grains, induc- 
ing shear banding in the model. Rollability of grains 
depends on: local shear, grain aspect ratios and local 
connectivity, which is related to the local density. There- 
fore, a straightforward extension to model shear banding 
and shear-rate-dependent boundary pressure would be 



by letting the local friction vary with the local density. 
Another straightforward way to include non-uniform dy- 
namic friction is by following phenomenological expres- 
sions suggested for ^id These extensions wiU be 
reported elsewhere. 

Note that, to model dense granular fluids, whose parti- 
cles are in contact all the time, we focused on quasi-static 
flows, where the stress relaxation in the fluid is faster 
than any other relevant physical process. Assuming in- 
stantaneous stress adjustment within plugs is particu- 
larly helpful for calculating the dynamics of coalescence 
of plugs, which has only been mentioned briefly here. It 
remains to be seen whether introduction of stress prop- 
agation dynamics (sound waves) is relevant to modelling 
dense granular flows as da Vinci Fluids. 
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